//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
/// \file SpinTracking/src/DetectorConstruction.cc
/// \brief Implementation of the SpinTracking::DetectorConstruction class

#include "DetectorConstruction.hh"

#include "G4RunManager.hh"
#include "G4NistManager.hh"
#include "G4Tubs.hh"
#include "G4Box.hh"
#include "G4Trd.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4VisAttributes.hh"

#include "G4UniformMagField.hh"
#include "G4FieldManager.hh"
#include "G4MagneticField.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4ClassicalRK4.hh"
#include "G4Mag_SpinEqRhs.hh"
#include "G4EqEMFieldWithSpin.hh"
#include "G4TransportationManager.hh"
#include "G4Transportation.hh"
#include "G4PropagatorInField.hh"

namespace SpinTracking
{

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4VPhysicalVolume* DetectorConstruction::Construct()
{
  // Get nist material manager
  G4NistManager* nist = G4NistManager::Instance();

  G4bool checkOverlaps = true;
  G4double cos_theta = 0.93;

  G4double BeamPipe_Rin =  3.15*cm , BeamPipe_Rout = 3.37*cm  , BeamPipe_z = BeamPipe_Rin/std::tan(std::acos(cos_theta))  *cm;   //BeamPipe_Rout = 3.37*cm
  G4Material* BeamPipe_mat = nist->FindOrBuildMaterial("G4_PLEXIGLASS");

  //
  // World
  G4double world_sizeXY = 2*BeamPipe_Rout;
  G4double world_sizeZ  = 2*BeamPipe_z;
  G4Material* world_mat = nist->FindOrBuildMaterial("G4_Galactic");   //

  auto solidWorld = new G4Box("World",                           // its name
    0.5 * world_sizeXY, 0.5 * world_sizeXY, 0.5 * world_sizeZ);  // its size

  auto logicWorld = new G4LogicalVolume(solidWorld,  // its solid
    world_mat,                                       // its material
    "World");                                        // its name

  auto physWorld = new G4PVPlacement(nullptr,  // no rotation
    G4ThreeVector(),                           // at (0,0,0)
    logicWorld,                                // its logical volume
    "World",                                   // its name
    nullptr,                                   // its mother  volume
    false,                                     // no boolean operation
    0,                                         // copy number
    checkOverlaps);                            // overlaps checking

  G4VisAttributes* worldVisAttr = new G4VisAttributes(G4Colour(1, 0, 0)); //
  worldVisAttr->SetVisibility(true);
  logicWorld->SetVisAttributes(worldVisAttr);


  // Beam Pipe
  G4ThreeVector pos = G4ThreeVector(0, 0, 0);
  auto solidBeamPipe = new G4Tubs("BeamPipe", BeamPipe_Rin, BeamPipe_Rout, BeamPipe_z/2, 0, 360 * deg);
  auto logicBeamPipe = new G4LogicalVolume(solidBeamPipe, BeamPipe_mat, "BeamPipe");
  new G4PVPlacement(nullptr, pos, logicBeamPipe, "BeamPipe", logicWorld, false, 0, checkOverlaps);

  G4VisAttributes* PipeVisAttr = new G4VisAttributes(G4Colour(1, 0, 0)); //
  PipeVisAttr->SetVisibility(true);
  logicBeamPipe->SetVisAttributes(PipeVisAttr);

  // Beam Pipe core
  auto solidBeamPipeCore = new G4Tubs("BeamPipeCore", 0*cm, BeamPipe_Rin, BeamPipe_z / 2, 0, 360 * deg);
  auto logicBeamPipeCore = new G4LogicalVolume(solidBeamPipeCore, world_mat, "BeamPipeCore");
  new G4PVPlacement(nullptr, pos, logicBeamPipeCore, "BeamPipeCore", logicWorld, false, 0, checkOverlaps);

  G4VisAttributes* PipeCoreVisAttr = new G4VisAttributes(G4Colour(1., 1., 1., 0.)); 
  PipeCoreVisAttr->SetVisibility(false); 
  logicBeamPipeCore->SetVisAttributes(PipeCoreVisAttr);

  // Set BeamPipe as scoring volume
  fScoringVolume = logicBeamPipe;

  //magnetic field
  G4double magneticFieldStrength = 1 * tesla; // 1T = 1.0 tesla
  G4ThreeVector fieldDirection = G4ThreeVector(0, 0 , -magneticFieldStrength); // -z direction
  G4MagneticField* BeamPipeMagFieldCore = new G4UniformMagField(fieldDirection);
  //spin tracking
  G4FieldManager* fieldManager = new G4FieldManager(BeamPipeMagFieldCore);
  logicBeamPipeCore->SetFieldManager(fieldManager, true); 

  //BMT eqution
  G4Mag_SpinEqRhs* pEquation = new G4Mag_SpinEqRhs(BeamPipeMagFieldCore);
    fieldManager->SetDetectorField(BeamPipeMagFieldCore);

    G4MagIntegratorStepper* stepper = new G4ClassicalRK4(pEquation, 12);

    G4double minStep = 0.01 * mm;
    // G4double minStep = 1 * mm;

    auto chordFinder = new G4ChordFinder(BeamPipeMagFieldCore, minStep, stepper);

    // Set accuracy parameters
    G4double deltaChord = 3.0 * mm;
    // G4double deltaChord = 10 * mm;
    chordFinder->SetDeltaChord(deltaChord);

    // Maximum allowed integration error in one integration sub-step
    G4double deltaOneStep = 0.01 * mm;
    // G4double deltaOneStep = 0.1 * mm;
    fieldManager->SetAccuraciesWithDeltaOneStep(deltaOneStep);

    G4double deltaIntersection = 0.1 * mm;
    fieldManager->SetDeltaIntersection(deltaIntersection);

    G4TransportationManager* transportManager = G4TransportationManager::GetTransportationManager();
    //------------------------------------------------addition by lvxr-----------------------------------------------------------------//
    G4Transportation* transportProcess = new G4Transportation();
    //transportProcess->EnableUseMagneticMoment(true);    // Old name - will be deprecated
    transportProcess->EnableMagneticMoment(true); // Whether to enable neutral particles to be deflected with force due to magnetic moment
    //--------------------------------------------------------------------------------------------------------------------------------//

    G4PropagatorInField* fieldPropagator = transportManager->GetPropagatorInField();

    // Limits for relative accuracy of integration
    G4double epsMin = 2.5e-7;
    G4double epsMax = 0.001;
    // G4double epsMin = 2.5e-5;
    // G4double epsMax = 0.01;

    fieldPropagator->SetMinimumEpsilonStep(epsMin);
    fieldPropagator->SetMaximumEpsilonStep(epsMax);

    fieldManager->SetChordFinder(chordFinder);


  return physWorld;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

}
